import numpy as np

def lu_decomposition(A, b):
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)

    n = len(b)

    # 初始化L和U
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    
    for i in range(n):
        # 计算U的第一行
        U[0][i] = A[0][i]
        # 计算L的第一列
        L[i][0] = A[i][0] / U[0][0]
    for i in range(1, n):
        # 计算U的第i行
        for j in range(i, n):
            U[i][j] = A[i][j] - np.dot(L[i, 0:i], U[0:i, j])
        # 计算L的第i列
        for j in range(i, n):
            L[j][i] = (A[j][i] - np.dot(L[j, 0:i], U[0:i, i])) / U[i][i]
    
    # 回代
    y = np.zeros((n, 1))
    for i in range(n):
        y[i] = (b[i] - np.dot(L[i, 0:i], y[0:i])) / L[i][i]
    x = np.zeros((n, 1))
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i + 1 : n], x[i + 1 : n])) / U[i][i]
    return L, U, x

# 测试提供的矩阵 A 和向量 b
A = [
    [30, 33, -43, -11, -38, -29, 37, 28, 23],
    [-480, -523, 644, 128, 621, 480, -618, -489, -329],
    [60, 266, -1862, -1991, 464, 546, -968, -1567, 1652],
    [540, 624, -782, 299, 493, 123, 567, 5, -122],
    [-450, -675, 2245, 2326, -1512, 1230, -822, 129, -189],
    [-300, -120, -1114, -1295, -1946, -302, -376, -1540, -609],
    [1080, 998, 508, 2460, 1628, 1358, 2896, 2828, 2002],
    [-1080, -1408, 3340, 2267, 21, -1202, 866, -2690, -1351],
    [-300, -435, 1594, 1685, 340, 2279, -27, 2917, -2336],
]
b = [[188], [-3145], [-4994], [680], [7845], [1876], [9712], [-11599], [10127]]

# 调用算法并打印结果
L, U, x = lu_decomposition(A, b)
np.set_printoptions(precision=2, suppress=True)
print("下三角矩阵L：")
print(L)
print("上三角矩阵U：")
print(U)
print("方程的解：", x.reshape(-1,1))

# 随机生成至少20阶的方阵和非零向量
n_min = 20
n = np.random.randint(n_min, n_min + 10)  # 生成20到30之间的随机整数
A_random = np.random.rand(n, n) * 20 - 10  # 随机生成的矩阵A
b_random = np.random.rand(n, 1) * 20 - 10  # 随机生成的向量b
L_random, U_random, x_random = lu_decomposition(A_random, b_random)
print("\n随机矩阵L：")
print(L_random)
print("随机矩阵U：")
print(U_random)
print("随机矩阵和向量求解的解：\n", x_random)

## 验证
# 验证解的准确性
b_hat = np.dot(A, x).reshape(-1)
b_orig = np.array(b).reshape(-1)
assert np.allclose(b_orig, b_hat, atol=1e-8), "解的验证失败：LUx ≠ b"

# 验证随机矩阵与向量的解
b_random_hat = np.dot(A_random, x_random).reshape(-1)
b_random_orig = b_random.reshape(-1)
assert np.allclose(b_random_orig, b_random_hat, atol=1e-8), "随机解的验证失败：LUx ≠ b"